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METHOD AND DEVICE FOR NUMERICAL ANALYSIS OF FLOW FIELD OF 
INCOMPRESSIBLE VISCOUS FLUID, DIRECTLY USING V-CAD DATA 

This is a National Phase Application in the United 

5 States of international Patent Application No. 
PCT/JP2003/015971, filed December 12, 2003, which claims 
priority on Japanese Patent Application No. 2002-379214, 
filed December 27, 2002. The entire disclosures of the 
above patent applications are hereby incorporated by 
10 reference . 

BACKGROUND OF THE INVENTION 
Field of the Invention 

The present invention relates to a method and device 
15 for numerical analysis of a flow field of incompressible 

viscous fluid. This numerical analysis directly uses V-CAD 
data that have substantial data including shape data and 
physical quantity data integrated into each other. 

20 Description of the Related Art 

Patent Literature 1 discloses a method for storing 
-substantial data. In this method, substantial data 
including shape data and physical quantity data that are 
integrated into each other can be stored with a small 
25 memory. Thereby, it is possible to manage the shape, the 
structure, the physical property and the history of an 



object in an integrated fashion. Further, it is possible 
to manage the data on a series of processes from designing 
to assembling, testing and evaluating in the same data so 
that CAD and a simulation can be integrated. 
[Patent Literature 1] 

Laid-Open Patent Publication No. 2002-230054 
As shown in FIG. 1, the method for storing 
substantial data including integrated shape data and 
physical property data of Patent Literature 1 includes an 
external data input step (Al) - £A-)-, an oct-tree division step 
(Bl) r and a cell data storing step (Cl) -fG-K At the 
external data input step (Al)- fAf , external data 12 
including boundary data of a target object that is obtained 
at an external data obtaining step SI are input to a 
computer or the like that stores the method. At the oct- 
tree division step (Bl) -fB-)-, the external data 12 is divided 
by oct-tree division into rectangular parallelepiped cells 
13 having boundary surfaces orthogonal to each other. At 
the cell data storing step (CI ) -£€-)-, various physical 
property values are stored for each cell. 

In the method of Patent Literature 1, the external 
data having shape data of the object is divided by oct-tree 
division into rectangular parallelepiped cells having 
boundary surfaces orthogonal to each other, and various 
physical property values are stores for each cell. The 
each divided cells are classified into an internal cell and 



a boundary cell. The internal cell is positioned inside or 
outside the target object, and the boundary cell includes 
the boundary. The internal cell has at least one kind of 
physical property value as attribute, and the boundary cell 
has at least two kinds of physical property values for the 
inside and the outside of the target object. 

The data treated in this method is called M V-CAD 
data", and a simulation using this data is called "Volume 
Simulation ,/ or "V-Simulat ion" . In FIG. 1, the reference 
numeral 14 designates V-CAD data. 

CFD (Computational Fluid Dynamics) has been gradually 
put into practical use. Accompanying this, grid generation 
requires more effort and time, and in the case of a 
complicated shape, the grid generation requires more time 
than the computation does. For this reason, fluid analysis 
using an orthogonal grid attracts attention. The fluid 
analysis using an orthogonal grid is described in Non- 
Patent Literatures 1 through 17. 

The experimental result on a flow around a forcibly 
oscillated circular cylinder" is described in Non-Patent 
Literatures 18. Calculation result by an ALE finite 
element method for self-excited oscillation caused by 
vortex generation from a circular cylinder is described in 
Non-Patent Literature 19. 
[Non-Patent Literature 1] 

Saiki, E.M. , Biringen, S . , 199 6, Numerical Simulation of 



a Cylinder in Uniform flow: Application of a Virtual 
Boundary Method, J. Comput . Phys .123,450-465. 
[Non-Patent Literature 2] 

Yabe Takashi et al, 1999 1997 , Solid-Liquid-Gas 
Unification Solving Method and CIP Method, Journal of Japan 
Society of Computational Fluid Dynamics, 7, 103-114. 
[Non-Patent Literature 3] 

Ye , T . , Mittal , R . , Udaykumar , H . S . , & Shvy, W.,1999, A 
Cartesian Grid Method for Viscous Incompressible Flows with 
Complex Immersed Boundaries, AIAA-99-3312 , 545-557 . 
[Non-Patent Literature 4] 

Nakano Akira, Shimomura Nobuo, Satobuka Nobuyuki, 
1995, Numerical Simulation of compressive Viscous Flows 
around an Arbitrary Shape Body Using Cartesian Grid System, 
Transactions of Japan Society of Mechanical Engineers, 61B- 
592, 4319-4326. 
[Non-Patent Literature 5] 

Ichikawa Osamu, Fujii Kozo, 2002, Computation of the 
Flow Field around Arbitrary Three-Dimensional Body Geometry 
Using Cartesian Grid, Transactions of Japan Society of 
Mechanical Engineers, 68B-669, 1329-1336. 
[Non-Patent Literature 6] 

PIAO Binghu, KURODA Shigeaki, 2000, Cartesian grid 
method for incompressible Viscous Fluid Flow, Journal of 
Japan Society of Fluid Mechanics, 19, 37-46. 
[Non-Patent Literature 7] 



Ono, K. , Tomita, N . , Fuj itani, K. , & Himeno, R. , 19 98, An 
Application of Voxel Modeling Approach to Prediction of 
Engine Cooling Flow, Society of Automotive Engineers of 
Japan, Spring Convention, No. 984 , 165-168 . 
[Non-Patent Literature 8] 

http : //kuwahara . isas . ac . jp /index . html 
[Non-Patent Literature 9] 

Teramoto Susumu, Fuji Kozo, 1998, Flow Simulation 
around Three-Dimensional Object Using a Cartesian Grid 
Method, Proceedings of 12 th Computational Fluid Dynamics 
Symposium, 299-300 . 
[Non-Patent Literature 10] 

Quirk, J. J. , 1994 , An Alternative to Unstructured Grids 
for Computing Gas Dynamic Flows Around Arbitrarily Complex 
Two-Dimensional Bodies, Computers Fluids, 23, 125-142. 
[Non-Patent Literature 11] 

Karman, S .L. Jr . , 1995, SPLITFLOW : A 3D Unstructured 
Cartesian/Prismatic Grid (12) ynamics of CFD Code for 
Complex Geome-tries, AIAA 95-0343. 
[Non-Patent Literature 12] 

Hirt,C.W.,& Nichols, B. D. , 1981, Volume of Fluid (VOF) 
Method for the D Free Boundaries, J. Comput . Phys . 39, 201-225 . 
[Non-Patent Literature 13] 

Hirt,C.W.,& Cook, J. L., 1972, Calculating Three- 
dimensional Flows Around Structures and Over Rough Terrain, 
J. Comput . Phys . 10, 324-340. 
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[Non-Patent Literature 14] 

Kase, Teshima, 2001, Volume CAD Development, Riken 
Symposium, Integrated Volume CAD System Research, The First 
Meeting, 6-11. 
[Non-Patent Literature 15] 

Toyoda, Arakawa, 1999, Analysis of Flow around 
Circular Cylinder Using Adaptive Cartesian Mesh Method, 
13 th Computational Fluid Dynamics Symposium, F03-1, CD-ROM, 
[Non-Patent Literature 16] 

Matsumiya, Koeda, Taniguchi, Kobayashi, 1993, 
Numerical Simulation of 2D Flow around a Circular Cylinder 
by Third-Order Upwind Finite Difference Method, 
Transactions of Japan Society of Mechanical Engineers, 59B- 
566, 2937-2943. 
[Non-Patent Literature 17] 

Bouard,R.,& Coutanceau, M . , 1 98 0 , The early stage of 
development of the wake behind an impulsively started 
cylinder for 40<Re<10 A 4, J. Fluid Mech., 101-3, 583-607. 
[Non-Patent Literature 18] 

Okamoto, S . , Uematsu, R. , and Taguwa, Y., Fluid force 
acting on two-dimensional circular cylinder in Lok-in 
phenomenon, JSME International Journal, B45, No . 4 , 
(2002) ,850-856. 
[Non-Patent Literature 19] 

Kondou, Numerical Simulation for Aerodynamic 
Behaviors of a Circular Cylinder, 15 th Computational Fluid 
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Dynamics Symposium, E09-2, (2001) , CD-ROM 

At present, in the fluid analysis, calculation of 
even a complicated flow field having a three-dimensional 
5 shape becomes possible by using an overlapped grid and 

unstructured grid method. However, mesh generation comes 
to occupy a large part of the entire simulation. For this 
reason, use of an orthogonal grid is desired for a mesh 
generation method that enables complete automatization. 

10 In numerical analysis on an arbitrary shape in an 

orthogonal grid system, it is generally difficult to treat 
an object boundary. Recently, several Cartesian grid 
methods have been proposed for discretizations near a fluid 
boundary, and a boundary condition. 

15 Specifically, there are proposed a virtual boundary 

of Non-Patent Literature 1, CIP (Cubic-Interpolated 
Propagation) of Non-Patent Literature 2, an immersed 
boundary method of Non-Patent Literature 3, NPLS 
(neighboring Point Local collection) of Non-Patent 

20 Literature 4, a method of introducing into a differential 
scheme a distance from the boundary located between the 
grid points of Non-Patent Literature 5, and a partial 
boundary adaptive Cartesian grid method of Non-Patent 
Literature 6. 

25 In these methods, the boundary of the object is 

strictly treated. However, to that extent, a computation 



process become more complicated. Accordingly, these 
methods are not necessarily suitable to a three-dimensional 
treatment for an arbitrary shape. 

On the other hand, in terms of practical use, two 
5 methods are promising. One method of the two forms a 
stepped boundary by using orthogonal grids so as to 
approximate an object shape (for example, Non-Patent 
Literature 7 of Ono in Nissan Automobile, Non-Patent 
Literature of Kuwabara in Computation Fluid Laboratory. 

10 The other method treats a boundary shape by introducing a 
cut cell to improve approximation (for example, Non-Patent 
Literature 9 of Fujii in Space Laboratory, Non-Patent 
Literature 10 of Quirk in J. J., NASA.). 

However, in the method using a cut cell, since the 

15 boundary extends through an arbitrary position in an 
orthogonal grid, cells neighboring each other on the 
boundary can greatly differ in size. For this reason, 
there is a report that viscous flow analysis is difficult 
in the cut cell orthogonal grid (refer to Non-Patent 

20 Literature 11) . 

When numerical analysis on a flow field of 
incompressible viscous fluid is performed by using the 
conventional overlapped grids and the unstructured grid, 
grid generation cannot be completely automatized. For this 

25 reason, the grid generation occupies a large part of entire 
simulation time, so that it was difficult to reduce 
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simulation time. 

Meanwhile, although numerical analysis on a flow 
field using an orthogonal grid enables the grid generation 
to be automatized, it is difficult to express the object 
5 boundary by using orthogonal cells. As a result, 

simulation accuracy becomes low. Particularly in the case 
of numerical calculation on a flow accompanied by the 
moving boundary, movement distance of the moving boundary 
is limited to integer multiples of a mesh that has a 

10 constant size, so that the calculation can become unstable. 

Furthermore, particularly in the cut cell method, the 
boundary extends through an arbitrary position in an 
orthogonal grid, so that the cells neighboring each other 
on the boundary can differ in size, so that viscous flow 

15 analysis was difficult in the cut cell orthogonal grid. 

SUMMARY OF THE INVENTION 
The present invention was made to solve the above 
problems. That is, it is an object of the present 

20 invention to provide a method and a device for numerical 
analysis of a flow field of incompressible viscous fluid, 
directly using V-CAD data, by which it is possible to 
completely automatize grid generation, to easily express an 
object boundary, and to perform simulation of high 

25 precision in a short time by a relatively simple 
calculation process without causing calculation 
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instability . 

According to the present invention there is provided 
a method for numerical analysis of a flow field of 
incompressible viscous fluid, directly using V-CAD data, 
5 comprising: 

a dividing step (A) of dividing external data into a 
plurality of cells (13) having boundaries orthogonal to 
each other, the external data including boundary data of an 
object which contacts incompressible viscous fluid; 
10 a cell classifying step (B) of classifying the 

divided cells into an internal cell (13a) positioned inside 
or outside the object and a boundary cell (13b) including 
the boundary data; 

a cut point determining step (C) of determining cut 
15 points in ridges of the boundary cell on the basis of the 
boundary data; 

a boundary face determining step (D) of determining a 
polygon connecting the cut points to be cell internal data 
for the boundary face; and 
20 a analyzing step (E) of applying a cut cell finite 

volume method combined with a VOF method to a boundary of a 
flow field to analyze the flow field. 

Further, according to the present invention, there is 
provided a device for numerical analysis of a flow field of 
25 incompressible viscous fluid, directly using V-CAD data, 
comprising : 
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an input device (2) for inputting external data 
including boundary data of an object (1) that contacts 
incompressible viscous fluid; 

an external storage device (3) for storing 
5 substantial data of shape data and physical property data 
integrated into each other, and a storage operational 
program for the substantial data; 

an internal storage device (4) and central processing 
device (5) for executing the storage operational program; 
10 and 

an output device (6) for outputting a result of the 
execution; 

wherein the device divides the external data into a 
plurality of cells (13) having boundaries orthogonal to 

15 each other, classifies the divided cells into an internal 
cell (13a) positioned inside or outside the object and a 
boundary cell (13b) including the boundary data, determines 
cut points in ridges of the boundary cell on the basis of 
the boundary data, determines a polygon connecting the cut 

20 points to be cell internal data for the boundary face, and 
applies a cut cell finite volume method combined with a VOF 
method to a boundary of a flow field to analyze the flow 
field. 

Furthermore, according to the present invention, 
25 there is provided a program for numerical analysis of a 

flow field of incompressible viscous fluid, directly using 
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V-CAD data, causing a computer to perform: 

a dividing step (A) of dividing external data into a 
plurality of cells (13) having boundaries orthogonal to 
each other, the external data including boundary data of an 
5 object which contacts incompressible viscous fluid; 

a cell classifying step (B) of classifying the 
divided cells into an internal cell (13a) positioned inside 
or outside the object and a boundary cell (13b) including 
the boundary data; 
10 a cut point determining step (C) of determining cut 

points in ridges of the boundary cell on the basis of the 
boundary data; 

a boundary face determining step (D) of determining a 
polygon connecting the cut points to be cell internal data 
15 for the boundary face; 

an analyzing step (E) of applying a cut cell finite 
volume method combined with a VOF method to a boundary of a 
flow field to analyze the flow field. 

By the method and device of the present invention, it 
20 becomes possible to perform simulation of high precision by 
a relatively simple calculation process, in which accuracy 
and stability can be obtained, calculation cost is not much 
required, a grid generation can be automatized completely, 
and an object boundary can be expressed easily. 
25 According to a preferred embodiment of the present 

invention, the analyzing step (E) comprises: 
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applying a two-dimensional QUICK interpolation scheme 

to a convection term for space integrals- 
applying central difference having precision of a 

degree of a second order to a diffusion term; 

combining the convection term and the diffusion term, 

and applying Adams-Bashf orth method having precision of a 

degree of a second order to the combined convection term 

and diffusion term for time marching; and 

applying an Euler implicit method having precision of 

a degree of a first order to a pressure gradient term for 

time marching. 

Thereby, it is possible to completely automatize grid 

generation with stability and discretization precision 

being maintained. 

For a two-dimensional boundary cell, a governing 

equation in the finite volume method is expressed by a 

governing equation (7), 

jj^idV = - jjdiv(u ®u)dV - jjdiv(pT)dV + — jjdiv(grad(u))dV ( 7 ) 



This equation (7) is obtained by integrating the 
equation (6) with a fluid part of the two-dimensional 
boundary cell being set as a control volume (CV) Vi,j. The 
equation (6) is a tensor divergence type rewritten from the 
basic governing equation (1) for the incompressible viscous 
fluid so that equation (7) can satisfy the basic governing 
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equation (1) for the incompressible viscous fluid. 

The convection term, the pressure gradient term and 
the diffusion term in the governing equation of the finite 
volume method are expressed by the equations (8), (9) and 
(10) , respectively, 

convection term: 

jjdiv(ii ®u)dV = (${u®u)ndS = £ -nSS m 

+*x{A iJ u% f2 ,, +I/2 v, . - , y -_ 1/2 v, )]/ ( 8 ) 

+[Av(H , v (x) u t j - B. , ,v (x) . ) 

-hAx(y4 / jVijV.j - v ij-* v i j-\ )]J \onfy no - slip on wall 
pressure gradient term: 

\\div{pI)dV = <#</>/) HdS= £ PJ'*SS m 

= Ay[B ijP ^ /2J -B^ jPi _ U2J ~P p (B,j -B^)]! ( 9 ) 

+Ax[4 J /?. y+1/2 -4.y-iA,/-i/2 -PpiAj -4. 7 -i)]7 
diffusion term: 

\\div{grad(u))dV = <^grad(u)-ndS = £ grad(u) m ^nSS ni 

= [Aj/(^ y gwirf(f/)* I/w -B^gradiiOl^. -(B.. -B^Jgradiii);) 
+M4jg*ri<MXj+v2 -Aj-ig^frXs-M-iAj - Aj^)g rad i ll Y P )^ (10) 

+Ax(4 ^ra^(v)^. +1/2 -A.j-yg^d(yXj-v2 -(4; -4,y-i)^^(v)^)]7 

The equations (8), (9), and (10) include object 
boundary data in which the boundary connects the cut points 
of the ridges so that it becomes possible to perform 
numerical analysis on the incompressible viscous fluid in 
the flow field at the boundary. 
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When a no-slip boundary condition is used for a solid 
boundary, integral is performed on the solid boundary with 
the convection term being zero, a value of a middle point P 
of a cut line segment being used as an average for the 
5 pressure gradient term and the diffusion term, and a space 
integral is performed with areas fractions being applied to 
all of the terms. 

Thereby, it is possible to maintain complete 
automatization of the Cartesian grid, and strictly satisfy 
10 the conservation law by the control volume in the fluid 
calculation . 

The boundary cell having the parameter smaller than a 
threshold value of VOF=0.01 is regarded as a complete solid, 
for the boundary cell having the parameter larger than the 

15 threshold value, a definition point for the parameter 

calculated in a cut cell is set at a center of the boundary 
cell, and a definition point for a parameter in a ridge is 
set at a center of a cell ridge, and a parameter at a 
middle point of a line segment 4 is calculated by a linear 

20 interpolation. 

Also in this manner, it is possible to maintain 
complete automatization of the Cartesian grid, and strictly 
satisfy the conservation law by the control volume in the 
fluid calculation . 

25 A drag force (in a flow direction) and a lift force 

(in a direction vertical to the flow) acting on the object 
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are expressed by equations (12) and (13), 
drag force : 

= JJ(%)^ + lli^dydx^^ds + ^ds (12) 



dx J J 3y 



^2 



only Cartesian 

lift force: 

v v °y s s 

y } x \ 

= l/ i( v) l/ 2 (v))^+ l ft( „ l fcW )*| 



only Cartesian 



These equations enable the drag force and the lift 
force to be calculated easily and accurately in the 
Cartesian grid. 

Preferably, in fluid-structure interaction analysis 
accompanying a moving boundary, a fluid system and a 
structure system are separately analyzed each predetermined 
time interval, and boundary conditions for the fluid system 
and the structure system are explicitly used. 

In this analyzing manner, the calculation process 
becomes simpler than that in the strong coupling that 
treats the fluid system and the structure system in an 
unified manner by coupling. Thereby, the calculation time 
can be reduced. Further, the movement distance is not 



17 



restricted to the integer multiples of the constant mesh 
size, so that avoid the calculation instability. 

In analysis on a forcibly vibrated circular cylinder, 
the circular cylinder is set as one-mass-point and one- 
5 degree-of-f reedom system such that the circular cylinder is 
a solid structure elastically supported and vibrating in a 
direction vertical to the flow, 

and Y-direction displacement of a center of the 
circular cylinder is given by the equation (17), and a 
10 velocity boundary condition in the Y direction for a 

surface of the circular cylinder is given by the equation 
(18) , 

y=Asin (2 iz f c t) ... (17) 

v w =A2 7U f c cos (2 7C f c t) ... (18) . 

15 Furthermore, movement velocity of the vibrating 

circular cylinder obtained by the equation (18) is changed 
to be given each calculation time step for the velocity 
boundary condition on the flow field. 

Thereby, it is possible to obtain the analysis result 

20 well consistent with the experimental data. For example, 
(1) the drag force and the lift force acting on the 
vibrating circular cylinder is larger than those acting on 
a stationary circular cylinder, (2) the lock-in phenomenon 
of the frequency of the vortex shedding due to the forcible 

25 vibration can be well reproduced at "0 . 15<f c <0 . 25" around 
the Strouhal number of 0.2, and (3) the drag force and the 
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lift force increase in the lock-in region. 

In analysis on self -induced vibration due to an 
vortex shedding from the circular cylinder, a vibration 
equation having a dimension is expressed by the equation 
(19) or (20), using one-mass-point and one-degree-of - 
freedom dumper /spring model, 

m^+c^ + ky = ±pU 2 0 DC L (19) 
+ (*r*/.)f&U (Ivfjy = £ C L (20) 

The movement velocity of the vibrating circular 
cylinder calculated by the equation (20) is changed to be 
given each calculation time step for the velocity boundary 
condition on the flow field. 

Initial displacement and initial velocity of the 
circular cylinder are set to be zero, the lift force is 
explicitly given by using a current value, and the 
vibration equation is integral by the Newmark' s (3 method 
to obtain vibration displacement and vibration velocity of 
the circular cylinder. 

In this analyzing manner, (1) the amplitude of the 
circular cylinder becomes the largest at the non- 
dimensional fluid velocity of Vr=4 (characteristic 
vibration frequency f o=l/4=0 . 25) , (2) the maximum 
displacement of y gradually increases in the range Vr=2 , 3 
and 4, and becomes the largest at the point of Vr=4, and 
(3) then, the maximum displacement of y decreases in a 
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waveform, accompanying increase in Vr . Thus, the result of 
this analyzing manner is qualitatively consistent with the 
calculation result of the ALE finite element method. 

Other object and advantageous features of the present 
5 invention will become apparent from the following 
description referring to the attached drawings. 

BRIEF DESCRIPTION OF DRAWINGS 
FIG. 1 is a flowchart showing a storage method for 
10 substantial data according to the prior application; 

FIG. 2 shows a configuration of a numerical analysis 
device for implementing a numerical analysis method 
according to the present invention; 

FIG. 3 is a flowchart of the numerical analysis 
15 method and a program therefor according to the present 
invention; 

FIGS. 4A and 4B are illustrations showing the 
definition of two-dimensional areas fractions; 

FIG. 5 is an illustration showing space 
20 discretization by a finite volume method at a boundary 
cell ; 

FIGS. 6A, 6B and 6C are images showing analyzed VOF 
distributions for three cases; 

FIGS. 7A, 7B and 7C are graphs for comparison on flow 
25 direction velocity distributions and theoretical solutions 
in three cases ; 
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FIG. 8 is an image showing an analysis grid and a VOF 
distribution; 

FIG. 9 is an image showing a velocity vector; 

FIG. 10 is an image showing a constant pressure line; 

FIG. 11 is a graph showing an average drag force 
coefficient according to the present invention; 

FIG. 12 is a graph showing average Strouhal number 
according to the present invention; 

FIG. 13A is a graph showing a drag force coefficient 
of a stationary circular cylinder according to the 
Cartesian adaptive grid method, and FIG. 13B is a graph 
showing average Strouhal number according to the Cartesian 
adaptive grid method; 

FIG. 14A is a graph showing a drag force coefficient 
of a stationary circular cylinder according to the general 
coordinate system boundary adaptive method, and FIG. 14B is 
a graph showing Strouhal number according to the general 
coordinate system boundary adaptive method; 

FIG. 15A is an image showing a flow in experiment at 
Re=300 and T=2.5, and FIG. 15B is an image showing a flow 
according to the present invention; 

FIG. 16A is an image showing a flow in experiment at 
Re=550 and T=2.5, and FIG. 16B is an image showing a flow 
according to the present invention; 

FIG. 17A is a graph showing temporal change in drag 
force and lift force according to the present invention, 
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and FIG. 17B is a graph showing the result according to the 
partial boundary adaptive Cartesian grid method; 

FIGS, 18A through 18F are images of application 
examples for practical use; 
5 FIGS. 19A through 191 show waveforms of temporal 

change in a drag force coefficient C D and a lift force 
coefficient D L acting on a circular cylinder vibrating at 
each frequency and a stationary circular cylinder; 

FIG. 20 shows relation between Strouhal number St 
10 expressing vortex shedding and forcible vibration 
frequency; 

FIGS. 21A through 21C show temporal change in 
displacement y (left drawing) , a drag force coefficient C D 
and a lift force coefficient D L (right drawing) , at a range 
15 of the non-dimensional fluid velocity of Vr=2, 3 and 4 
FIGS. 22D through 22G show temporal change in 
displacement y (left drawing) , a drag force coefficient C D 
and a lift force coefficient D L (right drawing) , at a range 
of the non-dimensional fluid velocity of Vr=5, 6, 7 and 8; 
20 and 

FIG. 23 shows the maximum displacement of y against 
the non-dimensional fluid velocity Vr . 
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DESCRIPTION OF THE PREFERRED EMBODIMENTS 
In the following, a preferred embodiment of the 
present invention will be described with reference to the 
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drawings . 

FIG . 2 is a configuration diagram of a numerical 
analysis device for implementing a numerical analysis 
method according to the present invention. The numerical 
5 analysis device of the present invention includes an input 
device 2, an external storage device 3, internal storage 
device 4, a central processing device 5, and an output 
device 6. 

The input device 2 is a keyboard, for example, and 

10 inputs external data 12 including shape data of an object 1. 
The external storage device 3 is a hard disk, a floppy disk, 
a magnetic tape, a compact disk or the like, stores 
substantial data constituted by a shape and a physical 
amount integrated to each other, and stores a storage 

15 operational program for the substantial data. The internal 
storage device 4 is a RAM or a ROM, for example, and stores 
operational information. The central processing device 
(CPU) 5 intensively processes the operation, the input, the 
output and the like, and executes the storage operational 

20 program in cooperation with the internal storage device 4 . 

The output device 6 includes a display device and a printer, 
for example, and outputs the stored substantial data and 
the executed result of the storage operational program. 
The numerical analysis device of the present 

25 invention employs the external storage device 3, the 

internal storage device 4, and central processing device 5 
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to divide the external data 12 into a plurality of cells 
having boundary surfaces orthogonal to each other, to 
classify the respective divided cells into an internal cell 
13a positioned inside or outside the object 1 and a 
5 boundary cell 13b including boundary data, to find a cut 
point in a ridge of the boundary cell 13b on the basis of 
the boundary data, to set as cell internal data the polygon 
that connects the found cut points, and to perform the 
analysis by applying to a flow field a VOF method combined 

10 with a cut cell finite volume method. 

FIG. 3 is a flow diagram showing the numerical 
analysis method and the program for this method. As shown 
in FIG. 3, the method and transformation program of the 
present invention includes a dividing step (A) , a cell 

15 classifying step (B) , a cut point determining step (C) , and 
a boundary face determining step (D) , and analyzing step 
(D) . 

The external data 12 input from the outside is: 
polygon data that expresses a polyhedron; a tetrahedron or 

20 hexahedron element used in the finite element method; 

curved surface data used in a three-dimensional CAD or a CG 
tool; or data that expresses surface of other three- 
dimensional shape by using information of partial planes 
and curved surfaces. 

25 Instead of these data (referred to as S-CAD data) , 

the external data 12 may be (1) data directly created by 
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manual input through an interface (V-CAD interface) 
specific to the V-CAD, (2) digitized surface data by a 
measurement device, a sensor and a digitizer, (3) volume 
data including internal information such as voxel data used 
5 in the CT scan, the MRI, and generally used in the volume 
rendering . 

At the dividing step (A) , the external data 12 is 
divided into a plurality of cells 13 having boundary 
surfaces orthogonal to each other. This external data 12 

10 includes boundary data of the object that contacts with 
incompressible viscous fluid, and is obtained at an 
external data obtaining step (not shown in the drawing) . 
In the case of a three dimension, this dividing is the oct- 
tree division, and in the case of a two dimension, this 

15 dividing is tetra division. 

In other words, at the dividing step (A) , a standard 
rectangular parallelepiped or a rectangle is divided into 8 
or 4 cells, and the dividing is recursively repeated until 
each cell is completely inside or outside the object. This 

20 dividing can greatly reduce a data amount compared to the 
voxel expression. 

One area formed by space division is referred to as 
the cell 13. The cell 13 is a rectangular parallelepiped 
or a rectangle having boundaries orthogonal to each other. 

25 A hierarchical structure, the number of divisions or 
resolution expresses an area occupied in the apace. 



25 



Thereby, the object in the space is expressed as the cells 
of different sizes that are put on top of each other. 

At the cell classifying step (B) , the respective 
divided cells are classified into the internal cell 13a and 
5 the boundary cell 13b. The internal cell 13a is positioned 
inside or outside the object, and the boundary cell 13b 
includes the boundary data. 

In other words, the cell positioned completely inside 
or outside the object that contacts with incompressible 
10 viscous fluid is classified into the internal cell 13a 

(cube) so as to have the maximum size. On the other hand, 
the cell including the boundary information in the external 
data 12 is classified into the boundary cell 13b. 

At the cut point determining step (C) , The cut point 
15 between the ridge of the boundary cell 13b and the boundary 
face is determined on the basis of the boundary data. 

At the boundary face determining step (D) , the 
polygon that connects the determined cut points are 
determined to be the cell internal data of the boundary 
20 surface. In the following, the cell that includes the 

polygon connecting the cut points is referred to as a cut 
cell. 

At the analyzing step (E) , the cut cell finite volume 
method combined with the VOF method is applied to the 
25 boundary on the flow field for the internal cell 13a and 

the boundary cell 13b. Thereby, the numerical analysis is 
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performed on the flow field of the incompressible viscous 
fluid. 

The present invention will be described in more 
detail . 

5 1. The present invention is aimed at a fluid 

analysis technique of practical use in a V-CAD program, and 
provides a method for analyzing an incompressible viscous 
flow around an arbitrary shape, on the basis of the cut 
cell (KTC) orthogonal grid. In the present invention, the 

10 boundary in the flow field is treated by using the cut cell 
finite volume method combined with the VOF (Volume of 
Fluid) method proposed by Hirt et al (refer to Non-Patent 
Literature 12). According to the method of the present 
invention, an internal flow that is the flow inside a 

15 channel, and an outer flow that is the flow around a 

circular cylinder are analyzed. The analyzed result was 
compared to the experimental data, the theoretical solution, 
and an analyzed result of the conventional method. 



20 2. Basic Equation and Calculation Method 

2.1 Basic Equation 

Basic governing equations are the Navier Stokes 
equation (1) and the continuity equation (2) shown in 
Formula 4 . 



25 



27 



[Formula 4] 



(1) 



dt dxj dx { Re dxfixj 




(2) 



5 



In Formula 4, 



Re" designates the non-dimensional 



Reynolds number defined by the characteristic length and 
the characteristic velocity of the flow field. This 
Reynolds number Re physically represents the ratio between 
the inertial force and the viscous force in the flow field. 

10 In Formula 4, "u" designates velocity, and "p" designates a 
pressure. "i = 1, 2, 3" and "j = 1, 2, 3" designate the 
respective directions in the orthogonal coordinate system, 
"i = 1, 2, 3" and "j = 1, 2, 3" are abbreviated to be "i" 
and "j" also in the following description. From the 

15 equation (1) , it is understood that the movement of the 
incompressible viscous fluid is controlled by only one 
parameter Re when an external force fi is not considered. 

2.2 Calculation Accuracy, Calculation Cost and 
20 Calculation Method 

Simulation on fluid is one kind of numerical 
experiment, and always accompanies a certain error. In 
order to accurately perform the fluid analysis, it is 
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necessary to satisfy the following 4 requirements: 

(1) Space resolution that becomes finer by acquiring 
the minimum length scale (such as a boundary layer, an 
vortex, a shock wave and a flame surface) of the flow. 
5 (2) A calculation area that becomes larger or the 

area where an artificial flow, influx/efflux boundary 
condition and blocking by the wall can be neglected, by 
sufficiently acquiring the maximum length scale of the flow. 

(3) Time and space discretization that enables a 

10 truncation error and numerical diffusion to be neglected. 

(4) Model suitable to a specific case (wall model, 
turbulent flow, and combustion model, for example) , and a 
scheme (K-K, QUICK, MUSCL, TVD, and ENO, for example) . If 
a finite difference method is to be used, the grid should 

15 be more orthogonal and uniform. 

The calculation cost (calculation time and necessary 
memory) of the fluid analysis can be determined by analysis 
accuracy required by the calculation method. The fluid 
analysis of practical use can be sufficiently satisfied by 

20 a certain accuracy (accuracy range of normal physics 
experiment, and accuracy required by a user) . 

Meanwhile, the flow field is classified into an 
internal flow, an external flow, and other flows (jet flow, 
for example) . Generally, an analysis area for the external 

25 flow is made larger than that for the internal flow. 

According to the present invention, in terms of practical 
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use, in the case of the external flow, all of the analysis 
areas are set to be 10D x 10D with the characteristic 
length scale of the flow field being D (=1) . 

According to the present invention, the finite volume 
5 method combined with the difference method is employed- In 
order to secure a certain level of stability and 
discretization accuracy, for space integral, two- 
dimensional QUICK interpolation scheme is applied to the 
convection term, and the center difference of the second 

10 order accuracy is applied to a diffusion term. For time 

advancement, the Adams-Bashforth method of the second order 
accuracy is applied to the convection term and the 
diffusion term, and the Euler implicit method of the first 
order accuracy is applied to the pressure gradient term. 

15 SOLA-HSMAC method (refer to Non-Patent Literature 13) that 
is proposed by Hirt et al is employed for coupling of the 
pressure and the velocity as an analysis algorism. In this 
case, repeated calculation is performed by using the SOR 
method with the easing coefficient being 1.65. The 

20 convergence is determined by a residual of 0.0002 in the 
continuity equation (2). Further, in order to prevent 
pressure vibration, defined points on the grid of three- 
directional velocities u, v, w, and a pressure p are 
positioned at a staggered grid that is displaced by the 

25 half mesh. 
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2.3 Treating of Solid-Fluid Boundary in V-CAD data. 
The most important algorism in the Volume-CAD system 
is the KTC algorism proposed by Kase and Teshima (refer to 
5 Non-Patent Literature 14), and restores the surface shape 
of the object from cut points of the cell. It should be 
noted that this KTC algorism has the same concept as that 
of the orthogonal grid cut cell method in the field of the 
fluid analysis. In the case of the KTC, even in the two- 

10 dimensional case, two or more intersection points (cut 

points) are present on 4 ridges of one cell. However, for 
simplicity, in the two-dimensional case, the number of cut 
points present on the ridges of one cell is limited to 2. 
In this case, the KTC has only 2 classifications as shown 

15 in FIGS. 4A and 4B. 

The areas fractions in the cut cell method are 
defined by the equation (3) of Formula 5, using the 
information of the intersection between the fluid-solid 
boundary and the cell ridge. 

20 In Foumula 5, "Axi,j" and M Ayi,j" designate the 

grid width in the X direction and the grid width in the Y 
direction, respectively, Ai,j and Bi, j designate the areas 
fractions in the X direction and the areas fractions in the 
Y direction, respectively, and Vxi, j and Vyi,j designate a 

25 line segment on the cell ridge occupied by the solid. 

By obtaining the areas fractions in the respective 
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direction, it is possible to obtain a volume void rate in 
the VOF method, using the equations (4) and (5) of Formula 
5. 

[Formula 5] 

Ax - — Vx. . Ay. . — Vv - 

Aj=— A — - : - °<4^, y .<i (3) 

VF U -1-0.5(2-4,, -A^d-B^-B^j) ( 4 ) 

when{A u j + A.j_ } > \)and{B u + B._ Xj > 1) 

VF Uj =0.5(4,, +*,-..,) ( 5 } 

w//en(4 ? , + 4,,^ < l)or(fi f - , + < 1) 



2.4 Discretization by the VOF finite volume method 

When conservative governing equation is discretized 
in the finite volume method, to apply the integral form of 
10 the equation (2) to the Green's theorem, the equation (1) 

in a differential form is rewritten into the diverging form 
of a tensor as shown in the equation (6) of Formula 6. 

In Formula 6, the mark of the combination of O and X 
designates a tensor product, and the vector I designates a 
15 unit tensor corresponding to Kronecker 5 i . 

The discretization at fluid-solid boundary is an 
important part of the present invention. For a grid other 
than a boundary cell, discretization in an orthogonal grid 
by the finite volume method does not cause a large problem 
20 even in a conventional technique. For this reason, the 
description of this discretization will be omitted. For 
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the fluid-solid boundary, assuming the KTC cell as shown in 
FIG. 4, space discretization is performed on the equation 
(6). The gradient of the velocity against time in a 
control volume is regarded as an average by using volume 
void rate so that time discretization can be performed on 
the equation (6) in a conventional manner. For this reason, 
the description of the time discretization in the boundary 
cell will be omitted. 

Assuming the fluid part of the two-dimensional 
boundary cell to be a control volume (CV) Vi,j, space 
integral is performed on the equation (6) (for simplicity, 
an external force term vector is x> f=0") . The governing 
equation in the finite volume method is expressed by the 
equation (7) of Formula 6. 
[Formula 6] 



The convection term, the pressure gradient term and 
the diffusion term in the equation (7) are applied to the 
Green's theorem (divergence theorem) so that the surface 
integral (for three dimension, volume integral) in the 
control volume is transformed into a line integral of a 
closed line (for three dimension, surface integral of a 
closed plane) . In other words, a line integral is 



riu 1 

~ = -div(u ® u ) - div(pl ) ^—-div{grad{u)) + / 



9t " Re 



(6) 




(7) 
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performed along the five ridges ( l—>2—>3-*4— *5—*l ) in the 
control volume shown in FIG. 5. For the Cartesian 
coordinate of the present invention, the convection term 
the pressure gradient term and the diffusion term are 
discretized as shown in Formula 7, respectively. 
[Formula 7] 
convection term: 

jjdiv(u<S>u)clV = <$(u<8>u)ndS = £ {S®u) m -BSS m 

V,.j V. m=1 " 5 

+Ax(4, . M ^> /2i , +1/2 v,, - Aj-^Uj-^j-v )V ( 8 ) 

^WX'U^- - b ^/1,^j ) 

+Ax(A iJ v£v iJ - A i } _ x v^v, j_, )]j\only no - slip on wall 
pressure gradient term: 

jjdiv{pI)dV = C$(pI)ndS = £ P J nSS m 

= Ay[Bj jP M , 2 J -^_ I ,,A-./2.y " " ( 9 ) 

diffusion term: 

ftdiv(grad(u))dV = <j$grad(u) ■ ndS = J] grad(u) m • fi5S m 

y>.i s.- 5 m= '- 5 

= W{B ijg rad(uy MI2j -B^gradWU,^ -(B^-B^gmdiu);) 
+Ax(A if grad(u)> j+V2 - A iJ _ l grad(u)> j _ U2 - (A t j - A f J _, )grad(u) p )]T (io) 
+iAy(B u grad(v)l U2J - B^gradiv)^ -(B tJ -B^grad^) 
+Ax{A i j grad(v)l^ l2 -A i j _,gradiy)l j _ V2 -(A, . - A iJ .,)grad{vY p )y 



It should be noted that the subscript of the veloci 
parameter varies depending on the staggered grid stencil. 
The vector i and j shows the equations in the X and Y 



34 

directions. The two-dimensional QUICK interpolation is 
applied to the velocities with the superscript (x) and (y) 
in the X and Y directions. The velocity gradients with the 
superscript x and y in the discretized equation (10) in 
5 Formula 7 designate the velocity gradients in X and Y 

directions, respectively, and the second order difference 
is applied to the respective velocity gradients with the 
superscript x and y. 

The discretization by the cut cell has two key points. 

10 One of the two key points is the manipulation of the 

integral in the solid boundary (cut line 4 in FIG. 5) . In 
the finite volume method, the boundary condition for the 
diffusion term of the second partial differentiation 
becomes the Neumann problem. In the case of no slip 

15 boundary condition in the solid boundary, the convection 

term is 0 at the cut line 4 in FIG. 5, so that the integral 
of the convection term is not necessary at the cut line 4, 
but the integrals of the pressure gradient term and the 
diffusion term are necessary at the cut line 4. In other 

20 words, the integral of the last term in each equation (9) 
and (10) is performed by using the middle point on the cut 
line 4 as the average value. It should be noted that for 
space integral, it is necessary to apply the areas 
fractions to the all terms. 

25 The other of the two key points is a defined point of 

a parameter calculated in the cut cell. In the finite 
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volume method, the volume average is obtained in the fluid 
void rate VOF, so that it is natural to position the 
defined point of the parameter at the center of the control 
volume of the fluid. However, if the defined point of the 
5 parameter is positioned at the center C indicated by the 
mark "O" having " • " therein of the control volume of the 
fluid, a complete automatization merit for the orthogonal 
uniform-space grid cannot be obtained. Therefore, 
according to the present invention, the threshold value of 

10 the VOF is set to be 0.01. In this setting, a boundary 
cell of which VOF is smaller than the threshold value is 
determined to be a complete solid, and for a boundary cell 
of which VOF is larger than the threshold value, the 
calculated parameter is not set at the control volume 

15 center, but set at the cell center (for example, the point 
D indicated by the mark "O" having " + " therein) . 
Similarly, the parameter defined points of the ridge lines 
3 and 5 in FIG. 5 are not set at the line segment center 
indicated by the mark "O", but set at the cell ridge line 

20 center indicated by the mark xx @" . Further, since a 

parameter is not defined at the center point P of the line 
segment 4, the parameter at the center point P is 
determined by the linear interpolation. In this manner, 
while the reproducibility of the solid boundary is 

25 deteriorated, it is possible to maintain complete 
automatization of the orthogonal grid generation. 
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Furthermore, in the fluid calculation, the conservation law 
can strictly be satisfied, and the above-described process 
is not likely to largely affect the calculation accuracy of 
the flow field. 

5 

2.5 Drag force and lift force acting on the object. 
At the time of analysis on the flow around the object, 
the calculation of the drag force and the lift force acting 
on the object becomes easy if the integral is performed 

10 along the first grid point of the object in the case of a 
boundary adaptive grid of the general coordinate system. 
Accordingly, in the present invention, the drag force and 
the lift force are calculated in the following manner. 

Closed integral is performed for fluid stress tensor 

15 acting on the surface of the object, and Green's theorem is 
applied so that the force acting on the object can be 
obtained as shown in the equation (11) in Formula 8, in the 
case of the two dimension. 

" o n " designates the stress tensor acting on the 

20 plane facing in the normal direction, and " a " designates 
the stress tensor at the surface of the object. It should 
be noted that the divergence of the second stress tensor in 
the two dimension is constituted by four directions. 
Accordingly, assuming the X direction as the flow direction, 

25 and the Y direction as the direction vertical to the flow, 
by the re-application of the Green"s theorem, the drag 
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force (in the flow direction) and the lift force (in the 
direction vertical to the flow) acting on the object can be 
calculated as shown in the equations (12) and (13) in 
Formula 8 . 

5 Once the drag force and the lift force are obtained, 

the drag force coefficient C D and the lift force 
coefficient C L can be made non-dimensional by using the 
main stream velocity U, the fluid density p , and the 
standard area A (in the case of the two dimension, the 
10 standard line) of the object, as shown in the equation (14) 
in Formula 8 . 

The non-dimensional Strouhal number St is calculated 
by the using the frequency f L of the C L as shown in the 
equation (15) in Formula 8. 
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[Formula 8] 

P 3<x rr 3<t . da 



S V V UJL ij V UJL j U Sj 



(11) 



drag force : 

= If(^)^ + JJ^T 1 )^ = j<r„ds+(j<r v eb (12) 
v ox v ay s s 

= UtF) UlJ->>*' + /(^ IftC) IfcW^lo* Cartesian 

lift force: 



CT „=-P + 2 A _ ^=a„=^- + -j + (14) 



^7 ^7 

Cd = P [/ 2 Z)/2 Ci = /7l/ 2 Z>/2 (15) 



St = f L DIU (16) 

3. Calculation result and consideration 
5 For testing an analysis scheme and a numerical value 

model in the numerical analysis on fluid, generally, the 
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benchmark test target for an internal flow is a channel 
flow, and the benchmark test target for an external flow is 
a flow around a circular cylinder. Also in the present 
invention, for testing the cut cell finite volume method 
5 combined with the VOF method, the numerical analysis was 
performed on the two-dimensional Poiseuille flow and the 
flow around the stationary circular cylinder, and the 
result of this numerical analysis is compared to the data 
obtained by other researchers. There is the analytical 

10 theorem solution for the two-dimensional Poiseuille flow, 
and there are much experimental data and calculation data 
for the flow around the stationary circular cylinder. For 
the application to the practical use, the analysis was 
performed on a flow over a backward-facing step, a flow in 

15 a narrow pipe, a flow around an obstacle in a channel, a 
flow in a branch pipe, a flow around a plurality of bluff 
bodies, and a flow in a premixed combustion chamber. The 
calculation result of this analysis is described in the 
following with reference to the visually shown drawings. 

20 

3.1 Testing on the channel Poiseuille flow 
For testing the effect of the VOF of the cut cell, 
analysis was performed on a duct in an analysis area of a 
10D x 10D square with the duct being inclined by 0 degree, 
25 10 degrees and 45 degrees. In this analysis, the number of 
the entire grids was 64 x 64, the half width of the duct 
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was 2.5D, the Reynolds number Re defined by this half width 
was 1, and the time interval At was 0.0001, the analysis 
time steps are 10000. The influx condition and the efflux 
condition both were set as a parabolic distribution. The 
5 pressure gradient in the flow direction was set to be 1.0. 
In addition, the HSMA method was used while directly giving 
the pressure boundary condition at the entrance, the exit 
and the wall. The VOF distributions for three cases are 
shown in FIGS. 6A, 6B and 6C. The broken lines in FIG. 6A, 

10 6B and 6C indicate the boundaries. The velocity 
distributions in the flow direction respectively 
corresponding to the three cases are shown FIG. 7A, 7B and 
7C, respectively with comparison to the theoretical 
analysis solution . 

15 In FIGS. 7A, 7B and 7C, the analysis indicated by the 

carved line A is theoretical analysis solution, and the 
"i=32" of the line curved line B is velocity distribution 
in the section vertical to the flow at the flow axis center 
positioned at the 32th grid point. The velocity 

20 distribution calculated for the duct that was not inclined 
(inclination angle: 0 degree) is well consistent with the 
theoretical analysis solution. Meanwhile, for the velocity 
distributions calculated for the ducts which are inclined 
by 10 degrees and 45 degrees, although the maximum velocity 

25 at the axis center is slightly underestimated, it can be 
said that the analysis result is consistent with the 
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theoretical analysis solution. 

3.2 Testing for the flow around the stationary 
circular cylinder . 
5 The analysis on the flow around the circular cylinder 

that is the testing target for the analysis scheme, the 
model would be the most difficult matter for the orthogonal 
grid. According to the present invention, the analysis on 
the stationary circular cylinder put in an uniform flow was 

10 performed for 24 cases in which the Reynolds Raynpldo number 
defined by the diameter of the cylindrical changes from 1 
to 20000, with the analysis area being 10D x 10D, and the 
grid number being 64D x 64D. The non-dimensional time in 
the calculation was set to be 100, and the coordinate of 

15 the circle center of the circular cylinder was set to be 
"X=2.5D, Y=5.0D" as shown in FIG. 6. For the velocity 
condition, the slip on the surface of the circular cylinder 
did not occur, the influx condition at the influx side was 
a uniform flow (velocity in the flow direction is 1.0), the 

20 efflux condition at the efflux side was a free efflux in 

the X direction (velocity gradient is 0) , and the free slip 
was set in the Y direction. Furthermore, since the 
velocity and the pressure are simultaneously adjusted by 
the HSMAC method, it is not necessary to directly give the 

25 pressure boundary condition for the influx, the efflux and 
the wall. 
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FIG. 8 shows the VOF distribution for the grid and 
the entire analysis area. FIGS. 9 and 10 show the velocity 
vector and the constant pressure line when the Reynolds 
number Re is 300. In FIGS. 9 and 10, the periodic Karman 
5 vortex street, the unsteady flow pattern, is observed 

behind the circular cylinder. This Karman vortex street is 
common when the Reynolds number Re is larger than 80, so 
that the calculation result for this in the case of the 
other Reynolds number is omitted. 

10 FIGS. 11 and 12 respectively show average values of 

the drag force coefficient and the average Strouhal number 
over the non-dimensional time 50 to 100, with compared to 
other experimental data. For the drag force, the 
calculation result for the Reynolds number Re smaller than 

15 60 in the present invention is well consistent with the 
result (B) obtained by the equation by Imai . Around the 
Reynolds number Re being 100, the calculation result of the 
present invention is slightly smaller than the experimental 
data (C) . For the Reynolds number being larger than 400, 

20 the calculation result of the present invention is not 
almost different from the experimental data (C) . 

The reason for this is considered to be that the 
resolution for the boundary layer at the wall at the 
circular cylinder is not sufficiently obtained when the 

25 Reynolds number Re becomes larger than the order of 100 
(Re=100(O)) because the space resolution in the present 
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invention is low (10/64=0.15). As understood from FIG. 11, 
for the Reynolds number Re being larger than 50, the 
calculation result (D) obtained by the step approximation 
without the cut cell is apparently overestimated to be 
larger than that of the experiment. On the other hand, the 
calculation in the cut cell VOF method proposed in the 
present invention greatly improves the overestimation for 
the Reynolds number Re being larger than 50. 

For the average Strouhal number, although the 
calculation result of the present invention is larger than 
the experimental data, the calculation result of the 
present invention stays within the range of 2.0 to 2.5 that 
is considered to be appropriate values for the Reynolds 
number Re larger than 200. 

FIGS. 13A and 13B show the analysis results of the 
average values of the drag force coefficient and the 
average Strouhal number for the circular cylinder according 
to the adaptive Cartesian grid method of other researchers 
(refer to Non-Patent Literature 15) . FIGS. 14A and 14B 
show the analysis results of the drag force coefficient and 
the Strouhal number according to the adaptive grid method 
for a general coordinate system (refer to Non-Patent 
Literature 16) . The Strouhal numbers obtained in these 
conventional methods are partially displaced from the 
appropriate range of 0.2 to 0.25 for the Reynolds number Re 
larger than 200. Thus, the calculation for the Strouhal 
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number according to the present invention (as shown in FIGS. 
11 and 12) is more precise than these conventional boundary 
adaptive grid methods. 

FIGS. 15A and 15B are enlarged illustrations of the 
5 flow around the circular cylinder showing the results of 
the visualized experiment of the Non-Patent Literature 17, 
and the calculation result of the present invention, for 
the Reynolds number Re of 300, at the non-dimensional time 
T of 2.5. FIGS. 16A and 16B are illustrations of the flow 
10 around the circular cylinder showing the result of the 
experiment of Non-Patent Literature 17, and the 
corresponding calculation result of the present invention, 
for the Reynolds number Re of 550, at the non-dimensional 
time T of 2.5. 

15 As understood from these drawings, in each of these 

case, the result of the present invention is well 
consistent with the visualized experiment for the shape and 
the size of the double vortex near the circular cylinder, 
and the distance from the directly back side of the 

20 circular cylinder to the position of the re-interflow of 
the circulated vortex. 

FIGS. 17A and 17B show temporal change in the drag 
force and the lift force of the circular cylinder according 
to the method of the present invention, and temporal change 

25 in the drag force and the lift force of the circular 

cylinder according to the c of Non-Patent Literature 6 by 
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PIAO et al. It was confirmed that the drag force and the 
Strouhal number obtained by the uniform interval Cartesian 
grid method of the present invention is well consistent 
with the result obtained by the partial boundary adaptive 
5 Cartesian grid method. The amplitude of the lift force 
obtained by the present invention is estimated to be 
smaller than that obtained by the partial boundary adaptive 
Cartesian grid method. It is considered that this is 
because the calculation area in the Non-Patent Literature 6 

10 is 30D x 15D, and a convection condition is used as the 
effluent boundary in the Non-Patent Literature 6. 

The calculation was performed for 24 cases of the 
Reynolds number Re of 1 to 20000 (not shown in the drawing) . 
As a result, the Karman vortex appeared at the Reynolds 

15 number equal to or larger than 80. This critical Reynolds 
number is somewhat larger than the critical Reynolds number 
Re of 50 in the experiment and the theory, but is 
considered to satisfy the requirement of the practical use. 

20 3.3 The application example for the practical use 

We attempted to apply the present developed analysis 
method to the several examples. As shown in FIGS. 18A 
through 18F, these examples are: (A) a flow over a 
backward-facing step; (B) a flow in a narrow pipe; (C) a 

25 flow around an obstacle in a channel; (D) a flow in a 
branched pipe; (E) a flow around a plurality of bluff 
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bodies; and (F) a flow in a premixed combustion chamber. 

In these drawings, only the u distribution for the X- 
direction velocities are shown. In these application 
examples, the calculation area was 10D x 10D, and the 
5 number of the grids was 64 x 64, similarly to the duct flow 
calculation (3.1) and the flow calculation around the 
circular cylinder (3.2). In these application examples, 
the Reynolds number Re was 100, and the influx condition 
and the efflux condition were set to be the same as in the 
10 calculation of (3.2). In other words, it becomes possible 
to calculate the flow field of an arbitrary shape by only 
inputting the V-CAD data directly, without amending the 
developed calculation program. 

15 4. Next, the analysis method of the present 

invention is applied to the numerical calculation (fluid- 
structure interaction analysis) on a flow accompanied by a 
moving boundary. The fluid-structure interaction analysis 
is important particularly from the engineering standpoint 

20 of the wind engineering and a heat exchanger. 

The numerical calculation was performed on (1) forced 
vibration of a two-dimensional circular cylinder placed in 
a uniform flow, and (2) a vortex self-excited vibration for 
a circular cylinder, as concrete examples of the fluid- 

25 structure interaction analysis. The resonance state caused 
by the lock-in phenomenon and the vortex self-excited 
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vibration was tested as a representative fluid-structure 
interaction problem that inflicts a fatal damage to a 
structure. The following description is directed to an 
effectiveness of the analysis method in the fluid-structure 
5 interaction analysis according to the present invention. 

4.1 Basic equation 

For calculation on circular cylinder vibration, an 
assumed model is a one-mass-point and one-degree-of -freedom 
10 system as an elastically supported rigid structure in which 
one mass point is vibrated in a direction vertical to a 
flow direction. 

In the case where the circular cylinder is forcibly 
vibrated, the circular cylinder is forcibly vibrated in a 
15 sine curve in a direction (Y direction) vertical to a 

uniform flow. In other words, the Y-direction displacement 
of the center of the circular cylinder is given by the 
equation (17) of Formula 9 where "A" designates an 
amplitude, "fc" designates the frequency of the forcible 
20 vibration, and "t" designates time. 

The velocity boundary condition in the Y direction on 
the surface of the circular cylinder is given by the 
equation (18) of Formula 9. 

In the case where the circular cylinder is forcibly 
25 vibrated due to the vortex shedding, the vibration equation 
having a dimension is given as the one-mass-point and one- 
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degree-of-f reedom dumper/spring model by the equation (19) 
of Formula 9. 

In the equation (19) , "m" designates the mass of the 
circular cylinder, "c" the viscous damping coefficient, "k' 
the elastic recovery coefficient of the spring, " p " the 
fluid density, "U 0 " the velocity of the uniform flow, "D" 
the diameter of the circular cylinder, and "C L " the lift 
force coefficient. The equation (19) is made non- 
dimensional by using the parameter, which controls 
vibration response, such as the characteristic frequency 
f 0 = ( k/m) °' 5 / (2 7z ) , the non-dimensional damping coefficient 
h=c/ (2 (mk) 0 . 5) , the Scruton number Sc=2 it h (2 p e / p ) , as 
shown in the equation (20) of Formula 9. 
[Formula 9] 

y = Asin(2jrf e t) ( 17 ) 
v w = A2xf e caB&rfj) 



m ^+ c ®+ k y = L pU l D C L 09) 

dt 2 dt * 2 y ° L 



^ + { 4xhf fl )^+ {Txffy - ^C t (20) 



4.2 Numerical value solution 

In the case where the circular cylinder is vibrated, 
the moving velocity of the circular cylinder calculated by 
the equation (18) or (20) is given as a value varying every 
calculation time step for the velocity boundary condition 
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of the flow field. 

For the analysis on the vortex self-excited vibration 
of the circular cylinder, the initial displacement and the 
initial velocity of the circular cylinder are set to be 0, 
and the vibration equation (20) is integrated by the 
Newmark' s j3 method to obtain the vibration displacement 
and the vibration velocity of the circular cylinder. 

The fluid-structure interaction analysis is not 
strong coupling that couples the fluid system and the 
structure system in a unified manner, but weak coupling 
that discretizes the fluid equation and the structure 
equation, and explicitly uses the boundary conditions for 
discretized equations. Particularly in the strong coupling, 
in principle, a time interval in coupling analysis on a 
complicated fluid-structure system is restricted to the 
smaller of the time interval for the fluid calculation and 
the time interval for the structure calculation. However, 
according to the present invention, the simple 
dumper/spring model is used in the calculation for the 
structure, so that the time interval is determined by the 
CFL condition in the fluid analysis. Further, when the 
future acceleration is integrated by the Newmark' s )S 
method, the lift force in the equation (20) is explicitly 
gived by using the current value. 
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4.3 The flow around the forcibly vibrated circular 
cylinder . 

In the case of the Reynolds number Re of 1000 in the 
flow field, the flow around the forcibly vaibrated circular 
5 cylinder was analyzed by using the equations (17) and (18), 
on the same condition as that of (3.2). It is known that 
for the Reynolds number Re of 1000, the Strouhal number St 
of the non-dimensional vortex shedding frequency for the 
stationary circular cylinder is about 0.2 as the 

10 experimental value. For the calculation condition on the 
forcible vibration, the amplitude A was set to be 0.1, and 
the calculation was performed for 8 cases, that is, the 
forcible vibration frequencies fc of 0.05, 0.1, 0.15, 0.2, 
0.25, 0.3, 0.4, and 0.5. FIGS. 19A through 191 show 

15 temporal varying waveform of the drag force coefficient C D 
and the lift force coefficient C L of the flow force acting 
on the circular cylinder at respective frequencies 
including the stationary circular cylinder. 

From the temporal waveform of the drag force and the 

20 lift force in FIGS. 19A through 191, it can be understood 
that the drag force and the lift force acting on the 
circular cylinder vibrating at the respective frequencies 
except for the frequency fc of 0.5 are larger than those 
acting on the stationary circular cylinder. 

25 FIG. 20 shows the relation between the forcible 

vibration frequency fc and the Strouhal number that 
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represents a vortex shedding state in the flow obtained by 
the vibration wave of the lift force. From the calculation 
result of fc/St=l at the range of 0.15<fc<0.25 around the 
Strouhal number of 0.2, it is understood that the lock-in 
phenomenon at the frequency of the vortex shedding caused 
by the forcible vibration is well reproduced. In addition, 
from FIGS. 19A through 191, it is understood that the drag 
force and the lift force increase in the lock-in region. 

The present invention and the experiment by Okamoto 
et al (of the Non-Patent Literature 18) differ (the 
Reynolds number (Re=19000) and the vibration amplitude 
(A=0.05) of the circular cylinder). However, the lock-in 
region of 0 . 186<f c<0 . 216 obtained by the experiment by 
Okamoto et at is included in the lock-in region of the 
present invention . 

4.4 The self-excited vibration due to the vortex 
shedding at the circular cylinder. 

On the same analysis condition as that of (4.3), the 
vibration displacement and the vibration velocity of the 
vortex self-excited vibration of the circular cylinder are 
calculated by integrating the vibration equation (20) by 
the Newmark' s j3 method. Further, the velocity boundary 
condition in the flow field is imposed for the vibration 
velocity of the circular cylinder. With the Reynolds 
number Re in the flow field being 10000, Scruton number Sc 
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being 5, and the damping coefficient h being 1%, the 
calculation was performed for 7 cases, that is, the non- 
dimensional flow velocities Vr=Uo/(foD) of 2, 3, 4, 5, 6, 7 
and 8. FIGS. 21A through 21C and 22D through 22G show the 
temporal changes in the drag force coefficient C D , the lift 
force coefficient C L and the displacement y for the 
respective the non-dimensional flow velocities. 

For the vortex self-excited vibration analysis on the 
circular cylinder shown in FIGS. 21A through 21C and 22D 
through 22G, the assumed model is one-mass-point and one- 
degree-of -freedom system of the elastically supported rigid 
structure that is vibrated in the direction vertical to the 
flow. However, from the calculation result shown in FIGS. 
21A through 21C and 22D through 22G, the lift force C L term 
at the right-hand side of the equation (20) irregularly 
varies due to the unsteady vortex shedding. Furthermore, 
for the fluid-structure interaction movement, the vibrating 
circular cylinder gives unsteady disturbance to the flow 
field. Similarly to the sine forcible vibration of (4.3), 
the averages of the drag force and the lift force caused by 
the vortex-induced vibration are larger than those in the 
stationary circular cylinder. Meanwhile, the time-varying 
curve is the waveform having larger irregularity than that 
in the sine forcible vibration of (4.3). The amplitude of 
the circular cylinder takes the largest value at the non- 
dimensional velocity Vr of 4 (the characteristic frequency 



53 



fo of the circular cylinder being 1/4=0.25). This is 
considered to be related to the Strouhal number St of 0.25. 

FIG. 23 shows the maximum displacement of y against 
each value of the non-dimensional flow velocity. The 
maximum displacement of y gradually increases at the range 
of Vr=2, 3, 4, and takes the largest value at the point of 
Vr=4 . Thereafter, accompanying the increase of Vr, the 
maximum displacement of y decreases in a waveform. It is 
confirmed that the varying tendency of the vortex-induced 
vibration response amplitude is qualitatively consistent 
with that in the calculation result by the ALE finite 
element method of the Non-Patent Literature 19 by Kondou et 
al. 

5. Conclusion 

(1) We developed the fluid analysis method that 
directly uses the V-CAD data, and performed the benchmark 
test on the internal flow and the external flow. As a 
result of this, we confirmed that the fluid analysis method 
of the present invention has the sufficient stability, is 
superior to the conventional voxel method for the 
calculation on the flow around the two-dimensional circular 
cylinder, and has the calculation precision that is not 
inferior to that obtained by the boundary adaptive grid 
method . 

(2) For the shape of the flow field, reading only 
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the V-CAD data enables the numerical simulation of the two- 
dimensional incompressible viscous flow to be performed on 
any complicated shape in a short time. 

(3) For the fluid-structure interaction vibration of 
5 the two-dimensional circular cylinder, in the case of the 
Reynolds number Re of 1000 in the flow around the circular 
cylinder, the lock-in region (fc/St=l) is expected to be 
the range of 0 . 15<f c<0 . 25 . 

Furthermore, it is demonstrated that the resonance 

10 state occurs due to the vortex self-excited vibration in 
the case of the Reynolds number Re of 1000, and the non- 
dimensional flow velocity Vr of 4, that is, when the 
characteristic frequency (fo=l/4) of the circular cylinder 
becomes equal to the Strouhal number (St=0.25) of the 

15 stationary circular cylinder. 

These calculation results are qualitatively 
consistent with the experimental data and the results 
obtained by the ALE finite element method. Accordingly, 
the fluid analysis method of the present invention is 

20 effective even in the fluid-structure interaction problem 
for the moving boundary. 
[Advantages of the Invention] 

As described above, we developed the numerical 
analysis method for the incompressible viscous flow field 

25 of an arbitrary shape in which the V-CAD voxel data 

expected to be the next generation CAD is directly used. 
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Furthermore, the cut cell finite volume method combined 
with the VOF method is used for the boundary of the flow 
field. Moreover, the method of the present invention was 
tested for the two-dimensional flows in the channel and 
around the circular cylinder. As a result of this, it was 
proved that the method of the present invention has a 
certain degree of precision and stability, and is a 
calculation method of practical use that does not spent 
much cost. 

Thus, the numerical analysis method and device of the 
present invention that directly uses the V-CAD data can 
achieve the excellent advantages in that it is possible to 
completely automatize the grid generation, to easily 
express the object boundary, and to perform the high- 
precision simulation in a short time by the relatively 
simple calculation process. 



